#ifndef PCPS_ACCUMULATOR_HEADER
#define PCPS_ACCUMULATOR_HEADER

#include <src/pcps.h>
#include <src/input.h>
#include <src/wavefunction.h>
#include <src/arena.h>
#include <src/coeff_value.h>

namespace pcps {

  // forward declaration of process accumulator
  template <class S, class PART, class REF> class ProcessAccum;

  //--------------------------------------------------------------------------------------------------
  // pcps::ThreadAccum -- Base class for objects that accumulate data during the evolution of one
  //                      thread's Markov chain.
  //
  //                      This class interacts closely with the ProcessAccum class, which
  //                      handles accumulation for a process.
  //
  //                      For a given flavor of accumulator (say that for the energy) each process
  //                      has one process accumulator and multiple thread accumulators.
  //                      The thread accumulators handle the actual accumulation of data during
  //                      the stochastic evolution of the threads' Markov chains.
  //                      The process accumulators provide a central storage area for the net
  //                      accumulations from all threads and communicate these accumulations
  //                      between processes.
  //
  //                      Thread accumulators are created from a process accumulator, and will be
  //                      linked to that process accumulator.
  //
  //                      Thread accumulators are derived from the Arena_Object class so that
  //                      their memory can be allocated using the Arena class.
  //                      See the Arena class description for why this is important.
  //--------------------------------------------------------------------------------------------------

  template <class S, class PART, class REF> class ThreadAccum : public pcps::Arena_Object {

    protected:

      // pointer to the corresponding process accumulator
      pcps::ProcessAccum<S, PART, REF> * const _pacc;

      // thread id
      const int _tid;

    public:

      // Declare, but do not define, a constructor.
      // The definition will come once the ProcessAccum class has been defined.
      // We use inline here to ensure the constructor has internal linkage.
      inline ThreadAccum(pcps::ProcessAccum<S, PART, REF> * pacc, const int tid);

      // destructor (declared virtual so derived classes' destructors are called when deleting parent class's pointer)
      virtual ~ThreadAccum() {}

      // function to reset the accumulator
      virtual void reset() = 0;

      // function to accumulate data during the markov chain evolution
      virtual void accumulate(const long int iter, 
                              const long int count,
                              const bool accept,
                              const bool finished,
                              const S coeff_ratio,
                              const int * const new_config,
                              const int * const old_config,
                              pcps::CoeffValueFunc<S, PART, REF> & cvf,
                              pcps::ThreadAccum<S, PART, REF> ** accumulators,
                              const int n_accumulators) = 0;

      // function to compute an update to the wavefunction (not all accumulators use this function)
      virtual void compute_update() {}

  };

  //--------------------------------------------------------------------------------------------------
  // pcps::ProcessAccum -- Base class for objects that accumulate results from all of a process's
  //                       Markov chains.
  //
  //                       For more info, see the description of ThreadAccum.
  //--------------------------------------------------------------------------------------------------

  template <class S, class PART, class REF> class ProcessAccum {

    protected:

      // data members
      const int _myrank;
      const int _nproc;
      const int _nthread;
      const long int _sample_length; // number of samples to take on EACH thread
      pthread_mutex_t _mutex;
      pthread_barrier_t _barrier;

    public:

      // data returning functions
      int myrank() const { return _myrank; }
      int nproc() const { return _nproc; }
      int nthread() const { return _nthread; }
      long int sample_length() const { return _sample_length; }
      double total_samples() const { return double(_nproc) * double(_nthread) * double(_sample_length); }
      pthread_mutex_t * mutex() { return &_mutex; }
      pthread_barrier_t * barrier() { return &_barrier; }

      // constructor
      ProcessAccum(const pcps::Input & userinp, const long int sample_length)
        : _myrank(MPI::COMM_WORLD.Get_rank()),
          _nproc(MPI::COMM_WORLD.Get_size()),
          _nthread(userinp.nthreads()),
          _sample_length(sample_length)
      {
        pthread_mutex_init(&_mutex, NULL);
        pthread_barrier_init(&_barrier, NULL, _nthread);
      }

      // destructor (declared virtual so derived classes' destructors are called when deleting parent class's pointer)
      virtual ~ProcessAccum() {
        pthread_mutex_destroy(&_mutex);
        pthread_barrier_destroy(&_barrier);
      }

      // function to create a thread accumulator linked to this process accumulator
      virtual pcps::ThreadAccum<S, PART, REF> * create_thread_accum(const int tid, pcps::Arena & arena) = 0;

      // function to reset the accumulator
      virtual void reset() = 0;

      // function to apply an update to the wavefunction (not all accumulators use this function)
      virtual void update_wfn(const pcps::Input & userinp, pcps::Wavefunction<S> & wfn) {}

      // function to print information during stochastic iterations
      virtual void print_iter_info(const pcps::Input & userinp,
                                   pcps::Wavefunction<S> & wfn,
                                   const int iter,
                                   std::vector<boost::shared_ptr<pcps::ProcessAccum<S, PART, REF> > > & accumulators) {};

  };

  //--------------------------------------------------------------------------------
  // pcps::ThreadAccum::ThreadAccum -- Now that the ProcessAccum class is defined,
  //                                   we can define this constructor.
  //--------------------------------------------------------------------------------

  template <class S, class PART, class REF>
  inline pcps::ThreadAccum<S, PART, REF>::ThreadAccum(pcps::ProcessAccum<S, PART, REF> * const pacc, const int tid)
    : _pacc(pacc),
      _tid(tid)
  {
    if (_tid < 0 || _tid >= _pacc->nthread())
      throw pcps::Exception( (boost::format("Error in ThreadAccum.  Thread id must be in the range [0,nthreads)")).str() );
  }

  //--------------------------------------------------------------------------------
  // pcps::OptSettings -- Structure to hold settings for optimization accumulators.
  //--------------------------------------------------------------------------------

  class OptSettings {

    public:

      bool opt_cor;     // whether we are optimizing the correlator variables
      bool opt_agp;     // whether we are optimizing the agp weight variables
      bool opt_orb;     // whether we are optimizing the orbital coefficient variables
      bool opt_par;     // whether we are optimizing the pairing matrix variables
      double step_size; // step size to scale the wavefunction update by
      int dim;          // dimension of the overlap matrix
      int n_var_corr;   // number of variable (i.e. non-fixed) correlators in the wavefunction
      int n_cor_var;    // number of correlator variables
      int n_agp_var;    // number of agp weight variables
      int n_orb_var;    // number of orbital coefficient variables
      int n_par_var;    // number of pairing matrix variables

      // constructor
      template <class S> OptSettings(const pcps::Input & userinp, const pcps::Wavefunction<S> & wfn) {

        opt_cor = userinp.optimize_corrs();
        opt_agp = userinp.optimize_agp_weights();
        opt_orb = userinp.optimize_orb_coeffs();
        opt_par = userinp.optimize_pairing_matrix();
        step_size = userinp.step_size();

        dim = 1; // starts at one because we include the old wavefunction as one of our basis vectors
        n_var_corr = 0;
        n_cor_var = 0;
        n_agp_var = 0;
        n_orb_var = 0;
        n_par_var = 0;
        if (opt_cor) {
          for (typename std::vector<pcps::Correlator<S> >::const_iterator c = wfn.correlators().begin(); c != wfn.correlators().end(); c++) {
            if (!c->fixed())
              n_var_corr++;
            if (!c->fixed() && !c->ti())
              n_cor_var += c->nelements();
          }
          dim += n_cor_var;
        }
        if (opt_agp) {
          dim += userinp.n_agp_weights();
          n_agp_var += userinp.n_agp_weights();
        }
        if (opt_orb) {
          dim += userinp.n_rest_orb_coeff();
          n_orb_var += userinp.n_rest_orb_coeff();
        }
        if (opt_par) {
          dim += userinp.n_pairing_variables();
          n_par_var += userinp.n_pairing_variables();
        }

      }

  };

} // end namespace pcps

#endif
